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Abstract 

We consider the problem of dispatching WindFarm (WF) power demand to individual 
Wind Turbines (WT) with the goal of minimizing mechanical stresses. We assume wind is 
strong enough to let each WTs to produce the required power and propose different closed-loop 
Model Predictive Control (MPC) dispatching algorithms. Similarly to existing approaches 
based on MPC, our methods do not require changes in WT hardware but only software 
changes in the SCADA system of the WF. However, differently from previous MPC schemes, 
we augment the model of a WT with an ARMA predictor of the wind turbulence, which 
reduces uncertainty in wind predictions over the MPC control horizon. This allows us to 
develop both stochastic and deterministic MPC algorithms. In order to compare different 
MPC schemes and demonstrate improvements with respect to classic open-loop schedulers, 
we performed simulations using the SimWindFarm toolbox for MatLab. We demonstrate that 
MPC controllers allow to achieve reduction of stresses even in the case of large installations 
such as the 100-WTs Thanet offshore WF. 
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1 Introduction 


In the last few years, the interest in wind energy has been constantly increasing. From the end 
2010 to mid-2013 the global wind capacity grew up by 48.3%, generating around 3.5% of the 
world electricity demand [I]. It has been estimated that, at the end of 2013, the worldwide 
wind capacity has reached 318 [GVF]. This constant increase of Wind Farms (WFs) installations 
is due to the fact that wind energy is an excellent environmental and friendly solution to the 
problem of energy shortage. For example, three years after the nuclear disaster of Fukushima, 
local Japan government, in particular the Fukushima prefecture, is considering to supply their 
regions with 100% renewable energy by 2040 [2]. To achieve this goal, despite the increase of the 
installed capacity of wind turbines, it is necessary to face new engineering and science challenges 
to improve efficiency and durability of the systems. In order to maximize the economic investment 
and the power generation efficiency, the size of WTs will be increased so as to produce more 
than 20 [MIF]. These larger WTs will be installed both in onshore and offshore environments, 
subjecting their flexible structures to forces of different entities. To face these problems, we need 
advanced control architectures: the aim is to improve the efficiency by reducing structural stress 
and hence extending the lifetime of components. Indeed, “the lifetime of wind power plants is 
considered to be about 30 years, even if usually after 20 years these plants are dismantled because 
of the progressive decrease in the energy production due to the aging of wind turbine components” 

i- 

WF control is essential to fit the required power, maximizing the performance, minimizing the 
mechanical forces acting on WTs and to detect anomalies in the WTs [3] and [S]. The required 
power is determined by a network operator, who specifies this value as a function of national 
load profiles and other economic and political criteria. However, the power that can be actually 
produced by a WF strictly depends on the wind blowing on its WTs. In this respect, each WT 
can work in two different operating regions. The first one is called power maximization region 
and it is selected when the wind acting on the WTs is not strong enough to ensure the production 
of the required power. The second one, called power tracking region, is selected if the wind 
blowing on the WT is enough to produce the demanded power. When all WTs are in the power 
tracking region, the use of a WF controller can bring major advantages, as one can choose different 
strategies to dispatch the power demand among WTs. A first simple solution is represented by 
the adoption of a scheduler: given the wind speed profile that acts on the WF and the power 
demand provided by the network operator, the scheduler divides the power demand according 
to an open-loop strategy, e.g. distributing the power equally between the WT or, alternatively, 
activating the smallest number of WTs. However, this could lead some WTs to be stressed much 
more than others. For these reasons, it is convenient to introduce a closed-loop WF controller, 
that uses the on-line measurements from the WTs (see Figure [IJ. 

In literature, there exist different approaches to the design of WF controllers. A first idea, 
explained in detail in [B] , is to use the knowledge of the available power Pa for each WT dispatching 
the FjXn proportionally to Pa- However, this approach could increase the tower oscillations and 
the stresses imposed on the motor shaft. In [7], the authors propose a linearized WT model and a 
Model Predictive Control (MPC) scheme that evaluates, on a given prediction horizon, the power 
demand set-points to achieve different aims such as mechanical stress reduction. An advantage 
of this controller is the possibility to force the fulfillment of input constraints. In this case, the 
stochasticity introduced by the wind is not considered as part of the model, but it is assumed 
that the wind profile is known a priori, leading to a deterministic approach. From [8], we can 
assume that the wind is the sum of an average speed, which change in the order of hours, and 
of a zero mean turbulence variation, that changes on a faster time scale. Using this assumption, 
low-pass filters are introduced in [siiiniin] in order to model the wind turbulence. In these papers, 
the authors propose stochastic WF controllers. However, they do not consider constraints on the 
references of power demand due to the mechanical characteristics of the WT and the instantaneous 
variations of the wind. It is worth noting that all the predictive control techniques discussed above 
allow to reduce the fatigue on the WF by making software changes on the SCADA systems and 
they do not require to replace hardware component of WTs. More recently, same authors of [7] 
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Figure 1: Wind farm controller. Inputs are the total power demand and, optionally, wind mea¬ 
surements. 


proposed a supervisory controller which can be easily installed on very large WFs [12] . For sake 
of completeness, it should be noted that in the literature there are examples of control schemes 
implemented directly on the WT (see for example m and M), but these approaches require 
substantial investments to upgrade existing WFs. 

In this paper, we propose new WF controllers using MFC regulators. As a reference model 
for a WT in the power tracking region we use a linearized version of the NREL model [15]. 
Differently from the previous approaches, we will account for the wind variations assuming a 
Kaimal wind turbulence spectrum [5] and modeling turbulence as an ARMA process. Then, we 
propose an optimal one-step-ahead predictor computed from the ARMA process. This allows us 
to develop Deterministic MFC (DMFC) and Stochastic MFC (SMFC) regulators with the goal 
of dispatching power demands between WTs so as to minimize tower bending and fatigue on the 
motor shaft while guaranteeing that the sum of the power demand for each WT meets 
particular, we will use the SMFC scheme proposed in [TB] and m in order to account for wind 
stochasticity and we also design two different DMFC regulators which will not account for the 
variance of wind turbulence. In this paper we will not make use of experimental data and, hence, 
to evaluate performance of the proposed MFC controllers, we perform several simulations using 
MatLab/Simulink and the SimWindFarm (SWF) toolbox [6]. 

The paper is organized as follows. In Section [5] we introduces the WF model, by proposing a 
linearized model of the adopted WT and an optimal one-step-ahead predictor for wind turbulence. 
In Section 131 we propose a SMFC regulator and two DMFC regulators. In Section Uj we present 
simulation results and Section[S]is dedicated to some conclusions and possible future improvements. 

Notation. We use a : b for the set of integers {a,a -I- 1,...,6}. The column vector with 
s components vi,...,Vs is v = (vi,... ,Vs). The function diag(Gi,..., Gg) denotes the block- 
diagonal matrix composed by s block Gi, i = l,...,s. Moreover, tr{Q) is the trace of matrix 
Q. The symbol 1^. denotes a column vector in R’’ with all elements equal to 1. Furthermore, I is 
the identity matrix. We use ||a;||p to define the P-weighted seminorm, defined for all x G M" by 
||a;||p = x^Px, where P is a positive-semidefinite real symmetric matrix. The functions E[-], uar[-] 
and cov[-] denote mean value, variance and covariance of random variables. The function std[v], 
where v = (vi, ..., Ug), denotes the sample standard deviation of measurements Vi, i = 1 : s. The 
function P(A) denotes the probability of the event A. The function WGN{a, /]) denotes White 
Gaussian Noise (WGN) with mean a and variance /3. The standard normal distribution with mean 
a and variance /3 is denoted with Af{a, /3). 
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2 WF model 


In this section, we propose a linearized model for the WF. We first introduce a linearized model 
of a WT operating in the tracking region and then we design an optimal predictor for the wind 
turbulence. 

2.1 NREL WT model 

The NREL WT model is an offshore 5-MW baseline variable speed wind turbine equipped with 
an active hydraulic pitch control. This model has been proposed in order to become a standard 
for large WTs. In this section, we derive a linearized model of the nonlinear system described in 
Figured In particular, we use the results described in [TB], taking advantage of simplifications 
introduced in [B]. We defer the interested reader to m and m for a complete description of the 
model. In the following sections, we describe each block of Figure [2j Moreover, a variable with 
index 0 means steady-state variable. 



Figure 2: Overview of the blocks of the NREL WT model. 


2.1.1 Aerodynamics 

The conversion of the wind energy in available green energy can be described by an aerodynamic 
model of the rotor, hence we can model how the energy captured by the rotor can be converted into 
driving torque of the rotating machine. In the NREL WT model this transformation is represented 
by the following static nonlinear equation 

Pa{t) = '^pR^vitfCp , with X{t) = (1) 

where Pa{t) is the wind turbine power [W], p is the air density [;^], R is the radius of wind 
turbine rotor [m], v{t) is the wind speed C'p is the power coefficient, X{t) is the tip speed 
ratio, /3(t) is the collective pitch angle [°], ujr{t) is the rotational speed of wind turbine’s rotor 
jradj parameter in (ED is a characteristic nonlinear function depending on tip speed ratio 

and pitch angle of the blade. Furthermore, defining the aerodynamic torque applied to the rotor 
shaft Tr{t) = and replacing Pait) using (P), we obtain 

Tr{t) = ^pRMtrCQiXit),m), CQ{X{t),P{t)) = (2) 
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where Cq is the torque coefficient. During the conversion process, part of wind energy is dissipated 
through a secondary effect that acts on the rotor of the wind turbine. This force operates per¬ 
pendicularly to respect to the rotor plane, producing a tower bending moment and, consequently, 
oscillations on the wind turbine. The force exerted is called thrust force {Ft{t)) and is modeled 
by the following nonlinear static relation 

m) = Mt{t) = hFt{t) (3) 

where Ct represent the thrust coefficient, h is the tower height, Mtit) is the tower bending moment 
caused mainly by the thrust force Ft{t), hence we do not consider any elastic force which could 
increase 

From and @, by linearization about the operating point vg, /3o, i^ro, Tro, Tfjo, we obtain the 
following linear models 

Tr{t) - TrO = KyTr{v{t) - Wq) + K^Tri^rit) - W^o) + KpTriPit) - j3g) (4) 

Mt{t) - Mto = KyMt {v{t) - Wo) + KojMt - Wro) + K^Mt {P{t) “ /^o)- (5) 

2.1.2 Wind turbine local controller 

Each WT is equipped with a local controller. The NREL WT local control system is simpler than 
other WT controllers: indeed, it does not use wind speed measurements and, moreover, does not 
provide additional blocks for oscillation damping (see for example [1] and references therein). The 
NREL WT control scheme consists of two tracking loops: the first to compute the power reference 
Prefit) and the second to compute the pitch angle reference prefit), based on the measure of 
the rotational speed Wg(t) of the generator. The NREL WT controller operates in the following 
configurations. 

• Power tracking. Prefit), boosted to compensate the generator efficiency and constraints on 
the generator rated power, tracks Pdemit)- Prefit) is set by a PI regulator where the error 
is computed as the difference of ojgit) to respect to the steady-state rotational speed of the 
generator Wgo- The nominal gains of the PI are adapted online based on Pdemit) and Prefit)- 

• Power maximization. Prefit) is fixed to zero and Prefit) is evaluated through a nonlinear 
function implemented in a look-up table. 

Furthermore, a switching logic alternates this two configurations under specific conditions. Since 
our aim is to control the set-point Pdemit) for each WT, in the following we will consider the 
power tracking configuration. In order to obtain a linearized model of the NREL WT, we need to 
study the static behavior of the WT. Since the WT operates in two different configurations, the 
static behavior is completely different. A detailed analysis is given in Section 4.2 of m- 
Linearizing the local adaptive PI controller described above around the operating point iPrefo, Prefo,(-^go), 
we obtain the following linear model of the PI regulator for computing Prefit) 

^iujfgit) - Ugo) = - Wgo) + ^(wg(t) -Wgo) (6) 

^(Prefit) - Prefo) = - ..go) - ^iu^git) - O^go) (7) 

at lyj 1 ijj 

where is a time constant of first order low-pass filter, that lumps the effects of the measurement 
device, a;7 is the filtered rotational speed of the generator and Kp and Kj are the gains of the PI 
controller corresponding to Prefo and Prefo- Note that, asymptotically, uj^ = ojg = Wgg. 
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2.1.3 Transmission 


The transmission system is a MIMO linear system describing the stiffness and damping of the low 
speed shaft generator, hence describing, through damped harmonic oscillators, the effects of 
and Tg (generator torque) on Ur, lOg and Ms (main shaft torque). This part of the system can be 
modeled as a shaft with lumped inertia, omitting the fast dynamics related to the shaft elasticity. 
Therefore, using results in [19], the following linear relations are obtained 


— {Ur{t) - COro) 
UJg{t) 

{Ms{t)-Mso) 


T I 2 7 Tro) ngb{Tg{t) Tgo)) 

rigbWrit) 

’^ghJr 




(Tgit) - Tgo) 


‘‘gb 


Jj- 'bl^uJ( 


■{Tr(t) — Tro) 


( 8 ) 

(9) 

( 10 ) 


where Ugb is the multiplication ratio of the gearbox and Jg and Jr are the inertia of generator and 
rotor, respectively. 


2.1.4 Pitch actuator 

The pitch actuator drives /3 to Pref- This variation is carried out via a servo drive that moves 
each blade on Pref- This set-point is reached using hydraulic pitch actuator. However, for the 
design of WF controller, we can assume that /3 ~ /3re/- 


2.1.5 Generator power controller 

In the NREL WT, the output electrical power P{t) is modeled by the static nonlinear equation 

P{t) = flUJg{t)Tg{t) 

where /i is the generator efficiency. Therefore the generator power controller can easily compute 
the generator torque reference Tg^^ as 


T^Ht) 


Prefjt) 

fJXVgit) ■ 


Moreover, a linearized model of equation m is 


T^^t) - Tgo = -^{Prefit) - Prefo) - ^{bOg{t) - UJgo). 


fiUJgO 




( 11 ) 


( 12 ) 


2.1.6 Electrical generator 

The generator dynamics is described by a lower pass filter. However, for the design of the WF 
controller, in the power tracking configuration, we can assume that Tg ^ Tg^f and Pout ~ Pref ~ 

Pdem • 


2.1.7 Linearized WT model 


input 

u^'^, the 


= iP, bOr, i 


— Tdem ^ 


= V-Vo 

WT 

y 

= {Mu M, 
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the linearized dynamics is given by 


^ ^ (13) 

yWT(^) ^ ^ + Df^d^^it) (14) 

where matrices and are obtained from (P | - (fTUl) and (fT^ . 


2.2 Optimal one-step-ahead predictor of wind turbulence 

The wind blowing on the wind farm generates an exogenous input v that acts differently on each 
WT: for this reason, during the design phase of the controller, it is important to use as much as 
possible the knowledge of the wind field (see also [20] and (l^). As common in the literature, we 
can rewrite the wind as f = h + h where v is an average speed, depending on weather conditions, 
which changes in the order of hours, and i; is a zero mean turbulence variation that varies on 
a faster time scale [HI P- 17]. This latter component of the wind is generally due to thermal 
conditions (e.g. variations in temperature) and the friction with the earth’s surface. A wind 
profile is characterized by the average speed and the turbulence intensity Tj defined as 

Ti = — 

V 

where Ui is the standard deviation. We can describe the turbulence spectrum as a function of 
frequency using the Kaimal spectrum m p. 23] and |5]), given by 


$5(w) = cr? 


(1 + 0 ;^) 

^ TTV 7 


where Ly is a length scale [m]. In the literature, it is common to simplify the wind model by 
assuming the wind speed variations are distributed as a WGN. In addition, each WT in a WF is 
affected by the presence of neighboring WTs. This effect is neglected in this study and will be 
considered in future research. 

In order to derive a linearized model of a WF, a linear model of the wind is needed. We identify 
an ARMA model for the turbulence variation and then obtain an optimal predictor. The ARMA 
process is described by a linear combination of previous outputs and previous inputs 

Wy{t) ^ WGN{Q,a'^) [3T]. We can identify and validate an ARMA process for wind profiles 
described by each pair {v, Tj). Let the ARMA process described by the transfer function G{z) = 
Then, the optimal one-step-ahead predictor m can be derived as 


V(z) 


C{z)-A{z) 

G{z) 


V{z) 


(15) 


where V(z) is the Z-transform of v(t) and V{z) is the Z-transform of the predicted turbulence 
variation v{t\t — 1). Minimal realization of (fT51) in the state-space yields to the model 

Xy{t -b 1) = AyXy{t) + ByV{t) 

V(t\t — 1) = GyXy{t). 

Moreover the prediction error is distributed as Wy{t). 


2.2.1 Importance of wind predictor 

In the following, we show through examples the advantages of using an optimal one-step-ahead 
predictor for the wind turbulence. We consider two set of measurements of wind speed: the first 
set is used to identify the ARMA process and the second one is used to validate the optimal 
one-step-ahead predictor The sets have been produced using the SWF toolbox which allows one 
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to generate wind profiles distributed according to the Kaimal spectrum. We identify ARMA 
processes by trying different combinations of previous measurements and Gaussian noise samples 
and we choose the optimal predictor that minimizes the Final Prediction Error (FPE). In the first 
example, we consider a wind speed described by u = 20 and T/ = 0.1, hence (t| = 4. We obtain 
the following optimal predictor for the wind speed 


Xy{t+ 1) 
v{t\t — 1) 


0 

0 

o' 



■f 

0.25 

0 

0 

Xy 

(t) + 

0 

0 

0.5 

0 



0 

0.9021 


0.4406 

0.5389] 


Ht) - 20) 


20 , 


(16) 


where estimated variance of the prediction error is 0.9010, i.e. it is distributed as lFGiV(0, 0.9010). 
The validation set and the predicted wind speed profiles are shown in Figure [S] 



Figure 3: Wind profiles for vq = 20 and Tj = 0.1: in blue the validation set and in red the 
predicted profile obtained using (HU). 


In the second example, we consider a wind speed described by u = 12 and Tj = 0.01, hence 
cr? = 0.0144. We obtain the following optimal predictor for the wind speed 


Xy(t + 1) 
v{t\t — 1) 


0.3613 

0.5 

[0.8885 


0.2621 

Q Xy[t) + 

—0.8208] Xy{t) 


1 

0 

12 , 


iv{t) - 12) 


(17) 


where estimated variance of the prediction error is 0.0036. The validation set and the predicted 
wind speed profiles are shown in Figure 0] In both examples, we note that variance of the original 



Figure 4: Wind profiles for vq = 12 and T/ = 0.01: in blue the validation set and in red the 
predicted profile obtained using (HZD. 
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wind turbulence v has been reduced by 75%. This is very useful in an MFC architecture, where 
we need to predict the behavior of each WT and hence to predict also the wind turbulence over a 
time horizon. 


2.3 WF model 

In this section, we derive a model of the WF. First, we derive a local model of a WT. We note 
that the dynamics in (fT51) and (1141) depend on the measurements of the wind turbulence. However, 
using MFC, we need to predict the wind turbulence. To this purpose, we will use the optimal 
predictor obtained in the previous section by setting v = vq where vg is the wind speed at the 
operating point (see Section [2.1.111 . 

Discretizing dynamics in (1131) and (fHl) with sampling time Tg = 1 sec0, augmenting the state of 
the WT using the states of the optimal predictor and using the fact that 

= V{t\t — 1) + Wy{t) = v{t), 


we obtain the following discrete-time LTI model 

+ 1 ) = + BaU^'^it) + BdaWyit) 

y 


where 


= CaX^'^it) -b DaV^'^it) + DdaWy(t) 


= (a 


7la = 


■JWT 

0 


B^Cy 

Ay -j- ByCy 


■ 

'bWT- 


rsyu 

, Ba = 

0 

5 Bda — 

d 

By 


c, = DY^Cy] , D, = Dda = Dr 


(18) 


and B^^, and D^^ discrete-time counterparts (obtained through 

exact discretization) of the corresponding matrices in (fTlTll and (ITTl) . 

In order, to derive a WF model consisting of N turbines, we need to group N WT models 
described by m- Therefore the WF mode0 is given by 


where 



x{t -b 1) = 

Ax{t) -b Bu{t) + Bdw{t) 



yit) = 

Cx{t) + Du{t) -b Ddw{t) 




u = 

. . . , ) 


y = iyr,---, 

vV). 

W = {Wy^l, 

. . . , Wy^]\[ ) 

w - 1FGV(0,S„), 

Eu, = diag(CT 

2 2 \ 
1, . . . , CTjy; 

A = diag(Ha,i, 


B = diag(Ha4, 

• ■ - ; Bd ^ ]\f ) 

Bd = diag(Hda.i 

; • ■ - ; Bda,N ): 

C = diag(Ga,i, 

■ • ■ , Ca,N) 

D = diag(Da,i, 

■ • ■ 1 jy), 

Dd = diag(DdQ4,.. 

■ , Dda,N)- 


(19) 


In the following section, we will use model (HU to predict the behavior of the WF. Moreover we 
will detail how we can take into account the wind measurements at each time-instant. 

^Note that the choice of the sampling time depends on the working frequency of the WF SC ADA system, that 
is usually IHz. 

^With abuse of notation, the state, input, output and disturbance variables of the 2-th WT, as well as matrices, 
are indicated with subscript i. 
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3 MPC regulators for WFs 

In this section we present different MPC regulators for achieving optimal power dispatching. We 
first introduce performance measures for assessing the quality of a dispatching algorithm. Then, 
we present the basic MPC formulation with chance constraints and, finally, we derive SMPC and 
DMPC regulators. 

3.1 Performance measures 

In order to evaluate performance of different regulators, quantitative criteria are needed. In this 
paper, we use the index J proposed in m 

J = Jp + Jms + Jivit (20) 


where 

• Jp is a measure of the power production and it is defined as 


Jp = 


\ 


N 


T/n NR 


rated 


{Piit) - Pref,i{t)f 


with Prated IS the wind turbine rated power; 

• Jms is a measure of the total main shaft fatigue and is defined as 


N 


Jm^ — 0.2 std 


i=l 


M,40 : T) 
2-106 


• Jmj is a measure of the fore-aft oscillation on the tower and is defined as 

N 


JMt = E 0.05 std 


i=l 


Mt40 : T) 
23 • 106 


( 21 ) 


( 22 ) 


Note that in Jm^ and JMt we use the standard deviation instead of using the rain-flow algorithm 
and Palmgren-Miner sum as in m- This is due to the fact that the working frequency of a WF 
SCADA system is not high enough to represent the damage fatigue obtained with a rain-flow 
count. This problem has been widely investigated in literature and we defer the interested readers 

to mm HI]- 


3.2 MPC formulation 

In the following, we use Xt, Ut, yt, Wt and dt instead of x{t), u{t), y{t), w(t) and d{t), respectively. 
At each time instant t, we solve the following MPC optimization problem over the prediction 
horizon Nh 




Ufc, \/k=t:t+Nh f— ' ^ 

k—t 

(23a) 

Xk+i = Axk -I- Buk + BdWk, yk = t :t + Nh, 

(23b) 

yk = Cxk + Duk + DdWk, yk = t :t + Nh, 

(23c) 

= 0, yk = t :t + Nh, 

(23d) 

Uk > <““) <p, ys=l:S,yk = t:t + Nh. 

(23e) 
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For short, in (1^ we used the index k (instead of the double index k, t) for referring to variables 
within the prediction horizon t : t + Nh- Note that cost function (I23all and constraint (1^ 
correspond to minimize the total load fatigue while guaranteeing that the total power generated 
by the WF fulfills the power demand required by the network operator. Moreover inequalities 
(I23cl) represent linear probabilistic input constraints, where Cg G K.^, > 0 (so that the input 

constraint is inactive for = 0), S' is the number of linear input constraints and p is a maximal 
probability of constraint violation. In view of (EU and (l 2 ^ . we set 


Q = diag 


QMt 0 

0 qMs 


<lMt 0 \ 

0 QMs\) 


with QMt 


0.05 0.2 

(23 • Nh “ (2- 106)2Ar^- 


Furthermore, we assume R = diag(ri,... ,rjv) and hence the only tunable parameters are > 
0, Vi = 1 : TV. 

In order to remove constraint (l23dD followine [5S1 d. 5371 Isee also fTU| and i), we look for a 
matrix T G that parameterizes the linear feasible set 

{uk G : iJfUk = 0} = {Tuk : Uk G 


This can be achieved using the following transformation in the input space 


Uk — T uj^ 


(24) 


where 


I 0 


0 


... O' 

0 

-I 1 


Differently from 0, m and next we show how to take into account wind measurements 
at time instant t. Since the value in (IT^ and ([13 is known for each WT, constraints (I23b|) 

and (I23cl) for fc = t can be rewritten as 


xt+i = AoXt + But + Bddt 
y{t) = CoXt + Dut + Dddt 


(25) 


where 


d = 




Aq = diag(Ao,i,..., ^o.w)) 


01 


0 A, ’ 


C = diag(C'o,i,..., Co, at) 

Co..= [C^^ 0]. 


Therefore, in (I23b|) the state Xt+i depends in a deterministic way on xt, Ut, since Wt is fixed. 
Summarizing, using (I23cp . (1^3 and (1^51) . we can rewrite the MFC optimization problem as 

t-\-N h 


min 

Uk, Vk^t-.t+Nh 

Xk+i = Axk + Buk 

yk = Cxk + Duk 

r{cJuk>uTn<p. 

where 

B = BT, D = DT, 

Note that constraint (I23bl) (resp. (I23cl) ') has 
(I26cll and (ESi)- 


T [ 2/fe Q + 1 Wfe y 

(26a) 

k—t 


Xt+i = AtjXt + Biit + Bddt, 

(26b) 

yt = CoXt + Dut + Dddt, 

(26c) 

BdWk, yk = t + 1 : t + Nh, 

(26d) 

DdWk, yk = t + 1 : t + Nh, 

(26e) 

Vs = 1 : 5, yk = t : t + Nh, 

(26f) 

R = T'^RT, cJ = cJt. 



been split into constraints (I26bl) and (I26dl) (resp. 
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3.3 SMPC regulator 

In this section we design an SMPC regulator. Our aim is to rewrite the stochastic problem (|26l) 
as a deterministic optimization problem solvable through Semi-Definite Programming (SDP) |26] . 
To this purpose we adopt the approach to SMPC proposed in [TB] and m- 

The optimization problem that must be solved online at each time instant t is 

min tr{MoPt)+ ^ tr{MPk) (27) 


with respect to the unknowns C/fe, Gk Pk, Uk, for k = t : t + Nh, Xk, ioi k = t + 1 : t + Nh, 9ks for 
k = t : t + Nh and s = 1 : S', Xk for fc = t -|- 3 : t -I- Nh and subject to the LMI constraints 


xt+i = Aoxt + But + Bddt 
Xk+i = Axk + Buk, yk = t + 1 : t + Nh 
Xt = Xt+i = 0 


Xt+2 = BdX-wBd ■ 


Xk+i 

AXk -t BGk 

B 

{*) 

Xk 


[ {*) 

{*) 



- 

Xk 0 


Xk 

Pk 

Gk 0 


Uk 


0 I 


0 


{*) 

i*) 


Xk 0 
0 E-1 

(.) 


> 0 , 


,>0 


^ 0 ? ^ — ^ + 3 :^ + Nh 


V/t — t + 2,..., Nh 


(28) 

(29) 

(30) 

(31) 


(32) 



'Xj 


Pj 

Uj 



0 


.{*) 

1 



> 0 


(33) 


Uk Gk 
Gl Xfc 
3 


>0, k = t : t + Nh 


-T " ^ max 9ks 

Ca Uk < —U.- 

“ 4 

9ks >0, k = t : t + Nh and s = 1 : S 

,T - . 1 1 

Cg UkC-s ^ Ok: 


k = t : t + Nh and s = 1 : S 


2 \erf~^{l — 2 p) 


k = t ■. t + Nh and s = 1 : S 


(34) 

(35) 

(36) 

(37) 


where {*) denotes the matrix transpose of the corresponding block in the upper triangular part, 
erf{-) is the Gauss error function and where 



Co^gco 

g’Sqd 

CSQDd 


'G'^QG 

C'^Qb 

G'^QDd' 

Mo = 

(.) 

T^RT + b'^QD 

b'^QDd 

, M = 

{*) 

T^RT + b'^QD 

b'^QDd 


L (“^) 

(*) 

Dd QDd_ 


[ ib 

ib 

P>I QDd_ 


The control law Ut is then obtained as 

ut = Ut- (38) 

A detailed derivation of problem (I27|) from problem (I26a|) is described in the Appendix. Here, 
we just highlight that, up to the linearization of a square-root function which is needed for getting 
the affine constraint ((35ll . feasibility of (l28)l - (l37)l implies that chance constraints (126111 are fulfilled. 
Moreover, the cost in (l?7l) provides an upper bound to the cost in (I26al) . As shown in [TB] . 
tightening of the constraints (12611) and relaxation of the cost in (I26al) are needed for recasting the 
original nonlinear optimization problem into an SDP problem. 
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3.4 DMPC regulators 

In order to design a DMPC regulator, we do not consider stochasticity in the optimization problem 
(Ell), hence Wk = 0,Vfc = t+ l : t + Nh, ■ Therefore, the MPC problem can be rewritten as 


t+Nh 


_ min y" IIVfcllQ + l|MfclU+l|efe.i:s||p 

Uk, Vk=t-.t+Nh ^ 

(39a) 

Xt+1 = AoXt + But Bddt, 

(39b) 

yt = Coxt + Ditt + Dddt, 

(39c) 

Xk+i = Axk + Buk, yk = t + 1 : t + Nh, 

(39d) 

Vk = Cxk + Duk, yk = t+l-.t + Nh, 

(39e) 

c^Uk < -1- ek,s, yk = t :t + Nh, Vs = I : S', 

(39f) 

£k,s >0, yk = t : t + Nh, Vs = 1 : S 

(39g) 


where the bar on a variable denotes the mean value. Moreover, we replace probabilistic constraints 
(l26fl) with linear constraints (I3M where we introduced the slack variables ek,s- Slack variables 
are also weighted in the cost function (IHHall . where we assume p > 0. Using a deterministic MPC 
regulator, we have to solve a QP problem at each time instant: from a computational point of 
view, even if the order of the LTI system (TT^ increases, the optimization problem can be solved 
online with high sampling rate |26| . Moreover, in absence of constraint on the inputs, constraint 
(I39fl) do not appear in the optimization problem (1391) . Hence, the optimal value of slack variables 
is €k,s =0, yk = t : t + Nfi, Vs = 1 : S', and we can solve (1!^ explicitly, obtaining 


= + 7^)-l {B'^QAxt + B’^QB^dt) (40) 


where 


Jit 

Ut+l 

Ut+Nh 



D 

0 

0 

0 

... O' 


Co 


Dd 


CB 

i) 

0 

0 

... 0 


CAo 


CBd 

B = 

CAB 

CB 

D 

0 

... 0 

, A = 

CAAo 

, Bd = 

CABd 


CA^i'-^B 

CA^’^-'^B 



CB D_ 






Q = diag(Q, 


TZ = diag(i?,..., i?). 


In the sequel, we will refer to this approach as Explicit DMPC (EDMPC). We note that in (I40|) 
the control inputs over the prediction horizon depend both on the measured state Xt and the 
wind turbulence measurements dt- Furthermore, (1401) can be easily implemented in a SCADA 
system without requiring optimization tools. On the other hand, since the input constraints are 
not involved in the MPC problem, the matrix R must be chosen properly, as we will show in the 
example section. 


3.5 On-line control actions 

Summarizing, at each time instant t, the power demand set-points for the WTs are computed as 

Pdejt) = Tu^P^it) + PdemO 


where Pdernii') — {Pdem,Pdem,N{J}^-i PdemO — (TdemO,!? PdemO,A^) ^'Ud, USing the 

receding horizon principle, [t) is the optimal value Ut obtained 
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• solving the SDP optimization (E71) for the SMPC regulator 

• solving the QP optimization (IH!B for the DMPC regulator 

• computing the control inputs (IMll for the EDMPC regulator. 

3.6 SWF controller 

In the simulation examples, we will compare performance of the proposed controller to the con¬ 
troller provided with the SWF toolbox (in the following SWFctrl). SWFctrl dispatches 
between the turbines proportionally to the available power at each turbine. In particular, the 
controller is based on the following equation 

Pa=Yl ( 41 ) 

where Pa,i is the available power, Vmeas,i is the measured wind speed and (7™“^ is the maximum 
power coefficient for the j-th WT. Therefore, the WF power demand is distributed as 

P . 

F^dem,i — F^dem p • 

Fa 

Note that dSD is the maximal value of Pa that can be obtained in O- 


4 Simulation examples 

In order to assess the performance of the proposed control schemes, we use the SWF toolbox 
[B]. The SWF toolbox allows to simulate a WF scenario using the Taylor’s frozen turbulence 
hypothesis. This hypothesis, illustrated more thoroughly in concerns the interactions among 
WT due to the wind. In all simulations we will not make use of simplifications introduced by 
Taylor’s hypothesis, therefore we do not introduce simplifications in generating an ambient wind 
field and we not reduce the complexity of wake effect models [6]. Moreover, a WF modeled 
using the SWF toolbox presents further nonlinearities that were not considered in the design 
of our controller, such as elastic forces for the tower bending moment and saturation for pitch 
actuator and WT local controller. In the following, we show the simulations performed using 
MatLab/Simulink. In order to solve online the SDP and QP problems, we have used YALMIP 
[28] and MOSFK [29]. 


4.1 WF composed of 10 WTs as in |Bj 


In the first example, we test our control architectures for a WF composed of 10 WTs arranged 
as shown in Figure [S] and proposed in [B]. For this example, we have set a WF power demand 
= 30 [MW], equally distributed by the scheduler on the 10 WTs, hence PdemO,i = 3 [MW], 
z = 1 : 10. Moreover, in the MPC cost function, we set = 0.06 for SMPC and DMPC and = 0.1 
for EDMPC and we require that \ui{t)\ < 0.1 [MWj. The parameters and the constraints on 
Ui{t) are set in order to guarantee good performance around the given set-point Pdem 0 ,i- Indeed, 
far from the set-point, predictions using a linearized model could be inaccurate. Moreover, the 
power demand set-points can be changed accordingly with the limitations given by the SCADA 
system. The prediction horizon is iV/j = 2. The wind speed at the operating point is vq = 12 [^] 
and its turbulence is Tj — 0.1. The state-space model of the wind optimal predictor, used for all 
WTs, is 


xT{t+^) = 
<{t) = 


0.7039 

0.5 

[0.4189 


-0.6178] 


"Co) 


(42) 


with variance of the prediction error equal to 0.3512. 

In Table[T] we summarize performance using different controllers. We note that using the SWFctrl 
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Figure 5: WF layout for example as in [6]. D1 = 600 [to], D2 = 500 [to] and D3 = 300 [to]. 
(Figure from [^). 


we achieve better performance in terms of tracking of the required power, however SWFctrl induces 
more mechanical stress, in particular for the main shaft. MFC schemes improve performance in 
terms of mechanical stress: indeed, compared with the open-loop controller, using MFC controllers 
we can improve performance at least of 19.69% for Ms and 3.29% fo Mt- We also highlight that 
performance of DMFC and SMFC are better than using EDMFC: this is due to the fact that the 
weights Vi are higher for EDMFC in order to guarantee that the power demand for each WT does 
not change more than 0.1 [MVF]. Ferformance of DMFC and SMFC are comparable: however, 
solving a QF has computational burden lower than solving an SDF. 



J 

Jp 

JMs 

JMt 

Scheduler only 

0.2551 

0.0027 

0.0611 

0.1913 

SWFctrl 

-140.77% 

10 . 64 % 

-574.71% 

-4.43% 

EDMFC 

7.26% 

7.44% 

19.69% 

3.29% 

DMFC 

8 . 96 % 

5.81% 

22 . 30 % 

4 . 75 % 

SMFC 

8.46% 

4.78% 

21.19% 

4.45% 


Table 1: Controllers performance for a WF composed of 10 WTs as in |^. Table entries have been 
obtained by averaging values obtained in 5 simulations of 15 minutes each. Top row: performances 
using open-loop scheduling. Other rows: percentage increment/decrement with respect to the 
values in the first row. Best performances are in bold. 

In FigureElwe compare performance of scheduler (open-loop strategy) and DMFC (closed-loop 
strategy) in a single simulation. We note that using DMFC we achieve two aims: i) we guarantee 
that the WF produces the power demand required by the network operator (FigureHJ by removing 
the power drop of the open-loop strategy; ii) we reduce mechanical stress, by reducing variations 
in Ms and Mt (Figures IBcl and [6d1 for the 9-th WT). We achieve our aims by changing the power 
demand set-points: in Figure [Sal (for the 9-th WT) we highlight that instead of a constant power 
set-point, we allow to change Pdem ,9 in a range of 0.1 [MID] that gives also good performance in 
the rate of change of the power demand set-point (usually < Q.l [MIF], see Figure ISbl). 

4.2 Performance using different prediction horizons 

In this section, we test the proposed MFC controllers in a WF composed of 3 WTs arranged as 
shown in Figured The conditions of the WF and the regulator parameters are equal to those used 
in Section [4. 11 With this example we aim at studying performance for different prediction horizons. 
Results are shown in Figure El We note that for all MFC regulators maximum performance are 
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(a) Power demand for WT 9. (b) Gradient of power demand for WT 9. 




(c) Main shaft for WT 9. (d) Bending moment for WT 9. 




(e) Wind on WT 9. 


(f) Output WF power. 


Figure 6: WF composed of 10 WTs: comparison between scheduler controller, i.e. open-loop 
strategy (red) and DMPC, i.e. closed-loop strategy (blue). 



D1 I D1 


Figure 7: WF layout for example using 3 WTs. D1 = 400 [m]. 


achieved with prediction horizons Nh = 2 and = 3. This means that, due to inaccurate wind 
predictions, performance decreases if the prediction horizon increases. Moreover we also note that 
designing an MFC controller with Nh — 0 corresponds to dispatching the power demand based 
on the knowledge of current state and current wind measurements only. Hence, the prediction is 
one-step ahead only. 
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(a) Improvements in terms of J. 


(b) Improvements in terms of Jp. 




(c) Improvements in terms of Jm^ • 


(d) Improvements in terms of Jmj • 


Figure 8: WF composed of 3 WTs: comparison between MFC controllers using different prediction 
horizons Nh- Each plot has been obtained by averaging results obtained in 5 simulations of 15 
minutes each. In all panels percentage of improvement with respect to the open-loop scheduler is 
shown. 


4.3 Performance without wind predictor 

In this section, we test the proposed MFC controllers in a WF composed of 3 WTs arranged as 
shown in Figure [T) We use WF conditions and regulators parameters as in Section O Moreover 
we set Nh = 3. In Figure |9] we show performance with and without using the optimal wind one- 
step-ahead predictor. We note that for J, Jp and Jm* the use of the wind predictor decreases 
the performance. However Jm^ increases, in particular using DMFC and SMFC. The reasons are 
the following: i) for the linearized output Mt^i we do not consider any elastic model of the tower 
oscillations that depend on the wind acting on the tower (see [18)1: ii) the optimal wind predictor 
is designed locally for each WT and hence it does not account for wind interactions (see (SO) 131]). 
These effects are more apparent if the wind turbulence increases. In future research we will also 
consider elastic model of tower oscillations and optimal wind predictors taking into account wind 
interactions among WTs. However, if our goal is to minimize main shaft fatigue only, the proposed 
wind predictors guarantee good performance. 

4.4 Thanet Offshore Wind Farm 

In this last example, we consider the Thanet Offshore Wind Farm [35], a WF composed of 100 
WTs arranged as shown in Figure 1101 For this example, we have imposed a WF power demand 
^dem — 300 [MW], equally distributed by the scheduler among the 100 WTs, hence PdemOi = 
3 [MW], i = 1 : 100. Moreover, in the MFC cost function, we set = 0.06 for DMFC and 
ri = 0.1 for EDMFC and we require that [^^(t)] < 0.1 [MW]. The prediction horizon is Nh = 3. 
The wind speed at the operating point is vq = 15 [^] and its turbulence is T/ = 0.1. For this 
example, we were not able to use SMFC since, at every time instant, it requires the solution to a 
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(a) Improvements in terms of J. 





(c) Improvements in terms of Jm^ 



^■Predictor j 
I INo Predictor] 


(b) Improvements in terms of Jp. 


^■Predictor 
I INo Predictor 


I 2- 

'■ I 

nl- 


(d) Improvements in terms of Jmj • 


Figure 9: WF composed of 3 WTs: comparison between MFC controllers using optimal wind one- 
step-ahead predictor (blue) and without using (red). Each barplot has been obtained by averaging 
results obtained in 5 simulations of 15 minutes each. In all panels, percentages of improvement 
with respect to the open-loop scheduler are shown. 



Figure 10: WF layout for the Thanet Offshore Wind Farm [35]. 


very large-scale SDP optimization problem. 

In Table [2] we summarize performance achieved by using different controllers. Compared with 
the open-loop controller and the SWFctrl, MFC controllers can diminish the mechanical stress of 
13.07% for Ms and 2.13% for Mt. 
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J 

Jp 

Jm, 

Jm, 

Scheduler only 

3.7429 

0.0024 

1.2708 

2.4698 

SWFctrl 

-0.34% 

0.01% 

-4.74% 

1.92% 

EDMPC 

4.07% 

0.15% 

9.18% 

1.45% 

DMPC 

5.84% 

0.22% 

13.07% 

2.13% 


Table 2: Controllers performance of the Thanet Offshore Wind Farm. Table entries have been 
obtained averaging 5 simulations of 15 minutes each. Top row: performances using open-loop 
scheduling. Other rows: percentage improvement with respect to the values in the first row. Best 
performances are in bold. 


5 Conclusions 

In this paper, we proposed MPC-based algorithms for dispatching a power demand for the whole 
WF among different WTs. The goal is to achieve minimization of the total mechanical stress. 
At the modeling level, we proposed to include in WT models a one-step ARMA predictor of the 
wind turbulence. We then demonstrated through simulations that this allows MFC dispatchers 
to achieve good performances in realistic scenarios. Future works will focus on increasing perfor¬ 
mances by improving the mechanical description of individual WTs as well as the model of the 
whole WF by accounting for interactions among WTs. 
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A Derivation of the SMPC problem (I271) - (p71) 


First, we recall the following results that will be useful in the sequel. 


Lemma 1 (Schur Complement). Let Z = 


' A 


B 

C 


he a symmetric matrix partitioned into blocks 


A, B, C where both A and C are symmetric and square. Assume that A is positive semi-definite 
and C is positive definite. Let S = A — BC~^B^ be the Schur Complement of C in Z. Then, 
Z > 0 if and only if S > 0. 


From (I26d|l . the mean value state dynamics can be obtained by neglecting Wk, and it is given 
by ([281), i.e. 

Xk+i = Axk + Buk, 'ik = t + l-.t + Nh, (43) 

where Xk = ]E[a;fc], iik = ]E[wfc]. Defining the error variable Sk = Xk — Xk and assuming a control 
law of the form 

Wfc = u/c + Kk5k (44) 

where Kk G one has that, for Vfc = t -I- 1 : t -I- Nh, 6k is zero-mean Gaussian random 

variable with covariance matrix Xk evolving as 


Xk+i = E{4+iCi} = (^ + BKk)Xk{A + BKkf + (45) 


Moreover, Xk ^ M{xk,Xk), ior Wk = t 1 ■. t Nh 1- 

Remark 1. We note that since it is possible to measure online fi, ojr, (through the sensors 

placed on each WT) and the states of the optimal predictor, we can assert that the state of the 
system at time t is measurable. Moreover, since we can also measure wind speed for each WT, we 
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can affirm that Xt+i is not affected by stochasticity (see also (I26bll ). Therefore, (14511 is initialized 
with (l29ll . i.e. 


Xi = Xt+1 = 0 

and, always according to (1051). we also have (lOOll. i.e. 

(46) 

Xt+2 = BdEyjBj. 

(47) 


We highlight that, since depends both from variables Xk and Kk, the dynamics of the 
covariance matrix is nonlinear. However, by relaxing constraint (145 1) from equality to inequality 
constraint and using Lemma [TJ we can rewrite (H5|) as dSID, i.e. 


Xk+i 

{*) 

[ {*) 


AXk + BGk 
Xk 
{*) 


Bd^ 


W 


0 


> 0 , 


k — t 3 t Nfi 


(48) 


where Gk = KkXk and (*) denotes the matrix transpose of the corresponding block in the upper 
triangular part. In other words, if there are Xk, Xk+i and Gk verifying (H51) . then one has 
Xk+i > 

Next, using (I26c|) and (|26ep . we rewrite cost function (|26all as 



Xt 


t+Nh 


Xk 


II 

Ut 

Mo 

+ ^ E 

II 

Uk 

Wlr 


_dt_ 




_Wk_ 



(49) 


where 



'G^QCo G^QD G^QDd 


'G'^QG 

C^QD 

C^QDd 

Mo = 

{*) T^RT + b^QD b^QDd 

, M = 

(*) 

T^RT + b^QD 

b^QDd 


. {*) {*) DlQDd_ 


L (*) 

(*) 

DjQDd 


Our next aim is to remove averages in the cost function (1491) . For this purpose, we proceed as 
described in m- We recall the following properties. 


E [X^FX] = E [tr (FXX^)] = tr (E [FXX ^]) = tr (FE [XX'^]) (50) 


E [XX^] = var [X] + E [X] E [X^] 

where X is a random vector and F is a square positive semi-definite matrix. 
Applying (1501) to (|00)l . we obtain 



( 


Xt 


Xt 

T~ 

tr 

MqE 


Ut 


Ut 



V 


dt. 


dt_ 



Nh 


ME 


Uk 

Wk 


Xk 

Uk 

Wk 


(•) 


(51) 
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Applying (ICTl to the highlighted part (•) we have 



Xk 


Xk 

T~ 


Uk 


Uk 



_Wk_ 


_Wk_ 



- 1 

?r 

_ 1 


1E[K 

Wfc wl]]+var 

1 - 

1 - 

1 _ 

Wk_ 




WfeJ _ 



E[a;fc] 


E[a;fc] 

1 ^ 

var[xk] 

cov[xk,Uk] 

cov[xk,'Wk] 

= 

E[Mfe] 


E[Mfe] 

+ 

cov[uk,Xk] 


var[uk] 

cov[uk,Wk] 


_E[wfc]_ 


_E[uifc] 


cov[wk,Xk] 

cov[wk,Uk] 

var[wk] 


E[a;fc] 


E[a;fe] 

T 

var[xk] 

var[xk]K'^ 

0 ■ 


= 

E[Mfe] 


E[uk] 

+ 

Kkvar[xk] 

Kkvar[xk]Kl 

0 



0 


0 


0 


0 

Eui 


_ 

1- 

H <S 

m N 


E[a;fe] 

E[ufc] 

T 

+ 

var[xk] 

Kkvar[xk] 

o' 

0 

uar[a;fc] 

0 

0 ■ ■ 
y;-i 

-1 r 

var 


0 


0 


0 

I 


J 



var[xk\Kl 

0 


0 

I 


Replacing E [xk], E [ttfe], var [xk] and Kkvar [xk] respectively with Xk, Uk, Xk and Gk, we obtain 



Xk 


Xk 

T' 


Xk 


Xk 

T 

'Xk 

O' 


'Xk 

0 ■ 

-1 


Gl 

o' 


Uk 


Uk 


— 

Uk 


Uk 

+ 

Gk 

0 


0 



0 

0 

I 


_Wk_ 


_Wk_ 



0 


0 


0 

I 



w 






(52) 


This relation is valid \/k = t + 2,... ,Nh, since, for these values of k, matrices Xk are positive- 
definite and therefore invertible. For the time instants j = t, t -I- 1, we have 



Xj 


Xj 

T' 


Xj 


Xj 


Uj 


Uj 


= 

Uj 


Uj 


_Wj_ 


w,_ 



0 


0 


(53) 


Let us now define 



Xk 


Xk 

T~ 


Uk 


Uk 



_Wk_ 


_Wk_ 



Then, relaxing the equality constraint (I52p . we obtain 


Xk 


Xk 

T 

1 

O < 


'Xk 0 ■ 

-1 

GT 01 

Uk 


Uk 

+ 

Uk u 


0 


=>5 

=l 

_0_ 


_ 0_ 


=1 

O 




L J 


Applying Lemma [1] we rewrite (ISTl) as the LMI (|5^ . i.e. 



'Afc O' 


Xk 

Pk 

Gk 0 


Uk 


0 I 


0 


Xk 0 



(*) 

0 E-1 

0 

(*) 

(*) 

1 


Vk = t + 2,...,Nh. 


(54) 


(55) 
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Similarly, for introducing Pj, j = t,t + 1, we obtain the LMI i.e. 



Xj 


Pj 

Uj 



0 


.(*) 

I 



(56) 


As a whole, an upper bound to the cost function in (I26all is provided by the cost in (l?7l) . 

Our last aim is to account for probabilistic input constraints (I26fll using the procedure proposed 
in m- In the following, for simplicity of notation, we neglect the index s and the time k appearing 
in (|26fD . Suppose we want to impose 


V{c^u > < p, 


(57) 


Note that 


r{c^u > u'"“^) = P 


C^U — E[c’^m] ^ 


a/ c^var[u]c \J c^var\u]c 


= l-V 


-E[c 


C" U — ItLIC^ul {iTnax 


< 


-E[c 


a/ c^var[u]c \J c^war[M]c 


Since u is given by (PI). where Sk is Gaussian, one has that the random variable 
distributed as A/*(0,1). Hence we can write 


-y/ c'^ var[u]c 


IS 


p{c^u > u™“^) = i-g 


{pnax _ 


y/WvaT\^^ I 

where g(x) is the standard Gaussian probability distribution. Therefore, we can rewrite (1571) as 


^max _ 

■\/ c^var[u]c 


> 1 — p. 


Note that G is strictly monotone and invertible. Hence C/ ^ is strictly monotone as well. Therefore 
we can state 

y/c^ var[u\c 

In conclusion we obtained the deterministic constraint 


E[c^?x] < c^var[u]^ Q ^ (I — p). 


This new deterministic constraint, which involves the expected value of the random variable c^u, 
can be rewritten using the Gauss’s error function erf{x), which verifies 

= i(l + er/(-^)) 

as 


Gonstraint 


E[c'^m] < (p"var[u\cj V2erf ^ (1 — 2p). (58) 

is equivalent to the existence of E[{t], var[u\ and ijj > 0 such that, simultaneously, 

e..rlu]c < (59) 

c'^E [u] < ^ (60) 
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Considering as optimization variable, instead of tjj, we note that it enters (EOl) in a nonlinear 
way. In order to obtain an affine constraint, we linearize about ~ ^ ^ (observe that 

0 < '0 ^ yiTiaa: j^gnce the linearization point lies in the middle of the interval). We get 




yTTiaX 

4 


0^ 

yTYiaX ' 


Summarizing, constraint (1581) is replaced with (1591) and 


(61) 


c^E \u] < 

4 


0^ 

yjmax ’ 


We highlight that since from dSl) the control input Uk depends on the Gaussian error Sk-, we have 
that Uk ^ WGN{uk,KkXk), for all k = t : t + Nh, and therefore the control law variance must 
be assumed as optimization variable. The expected value E)?!] is u, while the related variance 
depends on var[xk\ and Kk- In fact we have 


var[uk\ = var [uk + Kk6k] = var [uk + Kk{xk - E[xfe])] 
= var [uk + Kk{xk - Xk)] = var [Kk^Xk - Xk)] 
= Kkvar [{xk - Xk)] = KkXkKl 


and substituting Gk = KkXk we obtain 

var[uk] = KkXkX-^XkKl = GkX-^Gl. 
Now, in order to obtain an LMI constraint, we relax (EH) as 

var[uk] - GkX~^Gl > 0 
and applying Lemma [U we have the constraint 

”0“*' 01 >0, * = 

Ak_ 

Concluding, to manage probabilistic linear input constraints, we replace 


(62) 


with 


■dSll), i.e. 


■■t + Nh (63) 

k = t : t + Nh and s = 1 : S' (64) 

C UkCs < (—^TTi-Sw") ’ k = t:t + Nh and s = 1 : S (65) 

2 Ver/-i(l - 2p)J 

where IJk = var[uk], E [csUs] = CgUs and 9ks = '4’ks- 

Summarizing all above results, the optimization problem that must be solved online at each time 
instant t is (EZD-dSTl)- 


Uk 

iGl 


Gfc 

Xk 

3 


> 0 , 


C^Uk < -u; 


k = t 

^ks 
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